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Abstract 

Lifting operators play an important role in starting a lattice Boltzmann model from a given 
initial density. The density, a macroscopic variable, needs to be mapped to the distribution 
functions, mesoscopic variables, of the lattice Boltzmann model. Several methods proposed 
as lifting operators have been tested and discussed in the literature. The most famous meth- 
ods are an analytically found lifting operator, like the Chapman-Enskog expansion, and a 
numerical method, like the Constrained Runs algorithm, to arrive at an implicit expression 
for the unknown distribution functions with the help of the density. This paper proposes a 
lifting operator that alleviates several drawbacks of these existing methods. In particular, 
we focus on the computational expense and the analytical work that needs to be done. The 
proposed lifting operator, a numerical Chapman-Enskog expansion, obtains the coefficients of 
the Chapman-Enskog expansion numerically. Another important feature of the use of lifting 
operators is found in hybrid models. There the lattice Boltzmann model is spatially coupled 
with a model based on a more macroscopic description, for example an advection-diffusion- 
reaction equation. In one part of the domain, the lattice Boltzmann model is used, while in 
another part, the more macroscopic model. Such a hybrid coupling results in missing data at 
the interfaces between the different models. A lifting operator is then an important tool since 
the lattice Boltzmann model is typically described by more variables than a model based on 
a macroscopic partial differential equation. 

Keywords: Lifting operator, missing data, lattice Boltzmann models, macroscopic partial dif- 
ferential equations, hybrid models, Chapman-Enskog expansion, Constrained Runs, numerical 
Chapman-Enskog expansion. 

1 Introduction 

A lifting operator is, in a multiscale method, an important tool that maps macroscopic variables to 
microscopic/mesoscopic variables. In kinetic models, for example, a lifting operator will map low 
order moments, like the density p(x, t) that counts the number of particles in a point x e D x C K™, 
n e No to a distribution function f(x,v,t) that counts the number of particles in a point (x,v) 
in phase space where x e D x C K" and the velocity v e D v C K™. For practical applications one 
uses n € {1, 2, 3} and time t > 0. 

In these problems the macroscopic level is typically described by a few low order moments 
and their evolution is simulated by use of a macroscopic partial differential equation (PDE). For 
example, the evolution of the density p(x, t) can be represented by an advection-diffusion-reaction 
equation. While the microscopic/mesoscopic level is typically described by a Boltzmann equation 
that evolves the distribution function f(x,v,t). A lattice Boltzmann model (LBM) is a special 
discretization of the Boltzmann equation that is, for example, used to simulate complex fluid 
systems. 
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Some examples of complex flows for which LBMs are used are flows in complicated geometries, 
multiphase and turbulent flows. Applications can be found in [51 125) and a more recent review is 
given by Aidun et al. [I . The application of the lattice Boltzmann model to multiscale physics in 
fluids is discussed in |26j . Banda et al. [5] present a high order relaxation system for the multiscale 
lattice Boltzmann equation to obtain the incompressible Navier-Stokes limit. 

Kevrekidis et al. introduced a lifting operator to couple different scales in a dynamical system 
in the equation- free framework [17] . This allows to model the dynamics at the macroscopic level 
by using short bursts of the microscopic simulation. 

A problem of lattice Boltzmann methods is the determination of initial conditions, usually 
given by macroscopic variables. During initialization and spatial coupling in a hybrid model, a 
one-to-many map needs to be created, known as a lifting operator. In this article we discuss a 
hybrid LBM and PDE model that uses a lattice Boltzmann model in one part of the domain while 
another part of the domain is described by a macroscopic partial differential equation. These 
different levels of description create missing data at the interfaces that can be resolved with a 
lifting operator. 

Hybrid approaches have been formulated for various flow problems. A Lennard- Jones particle 
dynamics is coupled with a compressible Navier-Stokes in [12 . A Boltzmann, respectively a lattice 
Boltzmann, model is coupled to the Navier-Stokes equations in [TjJ] and [TS] . A LBM is also coupled 
with a Navier-Stokes in Peano, an adaptive mesh refinement framework with spacetree grids [20 . 
Furthermore, the Boltzmann equation is coupled to the Euler equations in [3]. 

Coupled models also play an important role in the simulation of materials. A review of 
atomistic-to-continuum coupling is found in [22 . A more detailed overview for coupling methods 
in hybrid models can be found in [13] and [11] , In [13] adaptive mesh and algorithm refinement 
is used in parts of the domain where a continuum description is replaced by a particle descrip- 
tion. Coupled molecular dynamics and lattice Boltzmann models based on Schwarz's alternating 
method is presented in [llj . Dimarco et al. [10] merge deterministic methods for the equilibrium 
part with particle methods for the nonequilibrium part and present results for the Boltzmann 
equation with Bhatnagar-Gross-Krook (BGK) approximation. 

The coupling — that will be discussed in this article — of LBMs with reaction-diffusion PDEs 
is earlier considered in [251 133] • 

We propose a general lifting operator that maps densities to distribution functions. It is illus- 
trated for a LBM, but we believe that it is applicable to general discretizations of the Boltzmann 
equation and it can map more moments to the corresponding distribution functions. 

The new method will be compared to the Chapman-Enskog expansion [6 , a well known an- 
alytical method for the initialization of a lattice Boltzmann model, and the Constrained Runs 
(CR) algorithm [30] . a numerical lifting operator. The Chapman-Enskog expansion writes the 
distribution functions as an analytical series of the density. The Constrained Runs algorithm is 
based on the attraction of the dynamics toward the slow manifold and expresses, in an implicit 
way, the unknown distribution functions with the help of the density in successive grid points. A 
numerical comparison of these methods is given in [33] for hybrid models that spatially couple a 
diffusion PDE model and a LBM. Although the Constrained Runs algorithm is very accurate, its 
major drawback is the computational expense. It achieves a high accuracy for the coupling of a 
lattice Boltzmann model with a diffusion- reaction PDE [29]. However, lifting in the CR-algorithm 
requires many additional LBM steps. This computational cost is too expensive to be useful in more 
complex problems in higher dimensions. These drawbacks became clear in [33] when comparing 
the different methods numerically. The intention of this paper is to alleviate them. 

In this paper we propose a numerical Chapman-Enskog expansion that seriously reduces the 
computational cost of the lifting. It combines the idea of the Constrained Runs algorithm with the 
Chapman-Enskog expansion and does not need an analytical derivation as the Chapman-Enskog 
expansion. This lifting operator is calculated before the simulation and finds the coefficients of 
the Chapman-Enskog expansion numerically. Once these coefficients are found the application 
of the lifting operator is just a stencil computation as cheap as the analytical Chapman-Enskog 
expansion. As a spin-off it also extracts the macroscopic PDE from the lattice Boltzmann model. 
This allows us to construct the hybrid model without deriving the macroscopic PDE analytically. 
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The numerical results show that the new lifting operator can also reach a high accuracy. Although, 
we illustrate and benchmark the new method on academic model problems, we believe that it is 
applicable to other discretizations of the Boltzmann equation. Furthermore, for the clarity of the 
presentation we have kept the boundary between the LBM and PDE domain fixed. In a real 
application this boundary might be moved adaptively, triggered by an error estimate similar as in 
adaptive mesh refinement. 

This work is organized as follows. In Section [2] the model problem is defined. It focuses, in 
particular, on a hybrid model that consists of a LBM in one part of the domain and a macroscopic 
equivalent PDE in another part. Section [3] gives an overview of existing lifting techniques that 
are used in the literature. The Chapman-Enskog expansion and the Constrained Runs algorithm 



are respectively considered in Sections |3.1| and 3.2 These methods are discussed in Section [373 



We tend to remove these drawbacks by considering a numerical Chapman-Enskog expansion in 
Section |4j Section [5] contains the numerical results. In Section [O] the proposed lifting operator is 
tested in a setting of restriction and lifting. The application of the lifting operator to the hybrid 



LBM and PDE model is considered in Sections 5.2 and 5.3 We conclude and give an outlook in 
Section [6l 



2 Model problem 

Kinetic models make use of the Boltzmann equation [25] that describes the evolution of a distribu- 
tion function f(x,v,t) (function space C^(D)) that counts the number of particles or individuals 
in point x E Dx C R™, n € No, with a velocity v € Dv C R", at time t > 0. The equation is 

-f(x,v,t)+v — f(x,v,t)+F(x,t) — f{x,v,t) = n. (1) 

This is an evolution law in phase space where F(x,t) is the external force and O an integral 
operator that models the reorganization of the velocity distribution due to collisions or other 
interactions. 

The collision operator can be approximated by a simpler Bhatnagar-Gross-Krook (BGK) model 
il = ui(f eq (x, v, t) — f(x, v, t)) [T5] in which the equilibrium distribution f eq (x, v, t) is given by the 
Maxwell-Boltzmann distribution |25j . The BGK approximation represents a relaxation towards 
equilibrium with an associated time scale r = l/u>. 

At the moment it is still computationally expensive to simulate or analyze a Boltzmann model 
numerically. The development of efficient numerical methods for models based on the Boltzmann 
equation is therefore an active research field. A possible increase in efficiency can be obtained by 
constructing a hybrid model. The kinetic model is then replaced with a macroscopic description 
in the regions of the spatial domain where this is justified, for example, away from reaction fronts. 
These macroscopic models are cheaper to simulate. 

Possible macroscopic models for fluid dynamics are described by the Navier-Stokes and Euler 
equations. One can derive both the Navier-Stokes as the Euler equations from the Boltzmann 
equation [31 ]. A lifting operator then transforms the variables of the PDE to the distribution 
function of the Boltzmann equation at the boundaries between the domains. 

In this paper we study a simple model that allows a detailed study of the lifting operator both 
for the initialization of the distribution function and in a hybrid context. The model uses a lattice 
Boltzmann model with an equilibrium distribution function that only depends on the density. 
While the macroscopic PDE is an advection-diffusion-reaction equation for the density only. This 
simple model allows a detailed analysis yet it is general enough to expect that the results can be 
extended to more realistic Boltzmann models. The correspondence of the lattice Boltzmann with 
the Boltzmann equation is discussed in (3H [21] • 

2.1 Lattice Boltzmann models 

A lattice Boltzmann model (LBM) [25J OH] is a special discretization of Eq. ([!]). It describes the 
evolution of one-particle distribution functions fi(x,t) = f(x,Vi,t) discretized in space x, time 
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t and velocity Uj. The velocities are taken from a discrete set defined by the geometry of the 
grid. The functions are represented as fi :^xT^l with X x T the space-time grid with 
space steps Axi in the direction of velocity Vi, time step At and T = {0, At, 2At, . . .} n [0, T]. 
Representation DdQq used for the description of LBMs stands for d dimensions and q velocity 
directions. D1Q3, for example, considers in a one-dimensional spatial domain only three values 
for the velocity vi = CiAx/At with Ci = z, i € { — 1,0, 1} the dimensionless grid velocities. 

The remaining of this section contains the description of the lattice Boltzmann equation in one 
dimension but can easily be generalized to more dimensions. 

The lattice Boltzmann equation (LBE) describing the evolution of the distribution functions 
(with BGK approximation and no external force in Eq. ([T])) is 

fi(x + CiAx,t + At) = {l-u)fi{x,t)+uf° q [x,t). (2) 

The equilibrium distributions are given by ft q {x,t) = \p{x,t), i e {-1,0,1} [27] in which the 
particle density p(x,t) is defined as the zeroth order moment of the distribution functions p(x, t) = 
o 1} fi( x 't)- These equilibrium distributions correspond to a local diffusive equilibrium. 
The focus of this paper concerns the initialization of a lattice Boltzmann model. Starting the 
LBM scheme from a given initial density includes some arbitrariness. The distribution functions 
fi(x 7 0) at time t = need to be constructed from a given density p(x, 0). When the initialization 
is not consistent it leads to solutions with steep initial layers |21) . 



2.2 Macroscopic models for LBMs 

This section contains descriptions that represent macroscopic equivalent PDEs specific for LBMs. 

Partial differential equations model, at a macroscopic scale, the evolution of the moments of the 
particle distribution functions like density p(x,t) = ^ 4 fi(x,t), momentum (f>(x,t) = '^2 i Vifi(x,t) 
or energy £(x, t) = \ J2i v fMx, t). 

The transition between the distribution functions and the moments is straightforward since 
the matrix M below is invcrtiblc. 




1 -1 | | /o I = M | /o | . (3) 



When we look at these functions in a point x at time t, they can be represented either as 
(/l)/o,/-i) T € M 3 or as (p,4>>0 T € ^ we focus on the complete discretization in space, 
with n spatial grid points, the function spaces are M 3xn . 

It can be shown that the diffusion PDE and the LBM are macroscopic equivalent [28] when 
considering D1Q3 and 

dp d 2 p 2 — luAx 2 „„, . 1 . 

£ =D d> D= ^AT' #W) = 3P(M), *€ {-1,0,1}. (4) 

This can be checked by using a Chapman-Enskog expansion. Here fi(x, t), i € {—1, 0, 1} is written 
as a series, each term containing higher order derivatives of p(x,t) [51 [2"81 15]. 

fi = fP+fMAx + fWAx* + f®Ax* + ..., (5) 

where 

h " h "3' h ~ 3wcV h ~ lSc^ 2 W" 
The macroscopic diffusion PDE Q is obtained by summing the series of the Chapman-Enskog 
expansion in ([5| over the velocities. This considers purely diffusive effects. 

For advection-diffusion problems with uniform velocity field a on D1Q3 the equilibrium distri- 
bution functions are given by [22] 

f! q (x 1 t) = -(l + — + -l *G {-1,0,1}, c 2 s = -c 2 , c=— , (6) 
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with the equivalent macroscopic description 

dp dp d 2 p 2-ujAx 2 

Similar results can be obtained for higher dimensional problems. In two spatial dimensions, 
represented by x and y, with equal space steps Ax = Ay, the macroscopic equivalence is given by 

dp „d 2 p „d 2 p _ 2 — lu Ax 2 „„„. . , 

i = D d£ +D W ^l^AF' f?(*,t)= WiP ( X ,t), t€{0,...,4}, 

*V ^0=*, t»i = ... = tfl 4 = ~, (8) 

for D2Q5 27 and 

^ +a — +a — = D— + D— D = 2 - - — 
<9i 3 cte w dy dx 2 dy 2 ' 3w Ai ' 

rear \ V i ' a (^i ' a ) 2 &2 \ r,i / \ 2 1 Ax 2 

fi 0M) = w iP\ 1 + —gr + 2?)' «e{0,...,8}, a=(o a; ,o v ), c s = g^j, 

/a;\ 4 1 1 

x=l 1, w = g, = -, ie{l,...,4}, w i=3g: i€{5, ...,8}, (9) 

for D2Q9 27]. 

With such analytical expressions from the Chapman-Enskog expansion available, a lifting op- 
erator can be constructed. Indeed, fi(x,t) is then written as a series in function of the given 
density p(x,t). This allows us to construct the distribution functions from a given initial density 
necessary to initialize the LBM. 

We will consider the lattice Boltzmann model as the 'exact' model and these PDEs of the 
density as the macroscopic approximations. From this, we can construct a hybrid model problem 
as outlined in Section [2~3l 



2.3 Hybrid models 

This section deals with the construction of the hybrid model problem. Consider the one-dimensional 
problem D1Q3 but bear in mind that a similar construction of a hybrid model can be done in higher 
dimensions. In particular, we couple a lattice Boltzmann model with the macroscopic equivalent 
PDE. The resulting hybrid domain for D1Q3 is presented in Figure [Tl A one-dimensional domain 
[a,b] is considered that couples the PDE Q on [a,l[ with the LBE ([2| on [l,b]. Furthermore, we 
assume periodic boundary conditions. Note that in both subdomains the same grid spacings are 
used in space (Ax) and in time (At) and the boundary I remains fixed for all times. Using a 
different space-time grid is of particular interest for future work but the same spacings are used 
to highlight the coupling error. Similarly, the boundary may be moved adaptively as in adaptive 
mesh refinement. In an actual physical problem it will be important to use the hybrid model that 
gives a Boltzmann description where shock waves, contact discontinuities or sharp gradients occur. 
Since these move in time it might be useful to work with a moving interface method as in [8 and 
[9]. However, this is not the focus of the current paper. 

On the domain [a, l[, we discretize the PDE Q with cell centered central differences in space 
and forward Euler time discretization. The grid points Xj with j £ {0, 1, . . .,£>} cover this domain 
and for these points it holds that 

p(xj,t + At) = p(xj,t) + (p(xj-i,t) - 2p(xj,t) + p(x j+ i,t)) . (10) 

For the grid points Xj with j & {p + 1, . . . , n — 1} in [/, b] the LBE ^ holds. 

The full domain has an initial condition p(xj, 0), Vj € {0, 1, . . . ,n — 1}. The lifting operator is 
required to formulate the initial conditions of the LBM domain fi(xj,0), Vj € {p + 1, . . . , n — 1}. 



5 







PDE domain 




LBM domain 










h 








o 






O 


O 


o 


o 








fo 








© 






© 


o 





o 






p(x ,t) p{x v ,t) 
















/-i 















O 


o 





o 


X-i = x n - 


-1 


X Xl x p 


Xp+l 






Xn-1 



a I b 

Figure 1: The domain [a, b] in the hybrid model is split into [a,l[ on which we solve the PDE 
model and [I, b] on which we solve the LBM. The solid points (•) represent the grid for the density 
p of the discrete PDE, the circles (o) represent the LBM variables (fi, fo, f-x) T ■ The periodic 
boundary conditions and the coupling are implemented with ghostcells which are drawn by dashed 
circles. The density in the ghostcells of the PDE domain, in x-i and x p +i, are found by taking 
J2i fi m %n-i and Xp+i, respectively. However, the ghostcells for the LBM domain, in x p and x n , 
require a lifting operator that lifts p to (/i, fo, f-i) T hi these points. 



The periodic boundary conditions lead to the following boundary conditions for the PDE 
domain: Vt : p{x^ u t) = £\ /,-(x n _i, t) and p(x p+1 ,t) = Y^i fi{ x p+i^)- 

The aim is to construct the boundary conditions of the LBM domain in such a way that Vi > 
and Vj € {0, 1, . . . , n — 1} the macroscopic density defined as 

/ ,\ fp&j't) if j'e{o,...,p}, 

p{Xj,t) = < (11) 

(Ej/ifei*) lf 3 e{p+l,...,n-l}, 

behaves as the density of a LBM solved on the full domain. 

To formulate these boundary conditions, a lifting operator is required that maps the density 
p(x,t) in the ghost points xq and x p , the unknown of the PDE, to the distribution functions 
fi(x, t), i G {—1, 0, 1} of the LBM. 

This can be generalized by considering a higher dimensional spatial domain. The remaining 
derivations in this paper focus on one dimension although they can also be generalized to more 
dimensions. 



3 Review of existing lifting operators 

This section gives an overview of existing lifting operators that map densities to distribution 
functions. An analytical expansion that expresses the distribution functions as a series of the 



density and its spatial derivatives is given in Section 3.1 while a numerical method is presented 
in Section [372] Section [373] discusses these methods which results in the motivation to propose a 
lifting operator based on the combined ideas of the analytical and numerical method. 



3.1 Chapman-Enskog expansion 



The Chapman-Enskog expansion [51 [2H1 IS], already discussed in Section 2.3 
lifting operator. It constructs a mapping from p(x,t) to fi(x,t), i € { — 1,0,1} (D1Q3). 



can be used as a 



3.2 Constrained Runs algorithm 

An alternative numerical procedure is the Constrained Runs algorithm discussed in this section. It 
is well known that in phase space the dynamics are quickly attracted toward a slow manifold |14j . 
For the problem we are studying the dynamics on the slow manifold can be parameterized by the 
density p(x, t). The distribution functions are then of the form {fx{p(x, t)), f (p(x, t)), f_x(p(x, t))}. 
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Because of Eq. ([3]), it is equivalent to determine {fi(p), fo(p), f-i(p)} or {p,(f>(p),£(p)}- The 
missing distribution functions {/i, /o, /-i} can be found by determining 4> and £ for a given p such 
that cj), £ and p lie on the slow manifold. This is the basic idea of the Constrained Runs (CR) 
algorithm that was proposed by Gear et al. [U] for stiff singularly perturbed ordinary differential 
equations (ODEs) to map macroscopic initial data to missing microscopic variables. It uses the 
numerical simulator to find the missing data such that the evolution is close to the slow manifold. 

This Constrained Runs algorithm can be applied to lattice Boltzmann models [3D]. The state 
of the lattice Boltzmann model can be split into 



u = (p) and v — 



where u € W 1 and v g ~R 2n for a LBM with n spatial grid points. The density p is known so u 
is given, while v is unknown since <p and £ are missing. Denote the known initial conditions as 
u(0) =u e E™. 

The idea is now to initialize v such that the evolution of v under the LBM is smooth of order 
to. The smoothness condition is defined by 



d m+1 v(t) 



= 0, 

t=o 



which is approximated by 



A m+1 v(t)^ A f"+ l d? " +1 y , (12) 

where A m is the well-known forward finite difference stencil on v(t), v(t + At), . . .. For m = 0, the 
converged v satisfies the smoothness condition, up to a certain tolerance, and it is an approximation 
to the point of intersection with the slow manifold. This is schematically represented in Figure 
[2] This iteration is always stable and the point of intersection is found to first order accuracy 
compared to the Chapman-Enskog expansion for the LBM with BGK collisions for one-dimensional 
reaction-diffusion problems [30]. For m > 1, multiple LBM steps are necessary to estimate the 
derivative. For to = 1, this is often interpreted as a backward linear extrapolation in time [32j . 

The number of LBM steps used in the backward extrapolation determines the accuracy of the 
scheme. Higher order schemes increase the accuracy but they can become unstable. In this 
instability is circumvented by formulating the point of intersection as a fixed point 

v k+1 =C rn (u Q ,v k ), (13) 

where C m denotes one step of the CR-algorithm and to is related to the order of the time derivative 



that is set to zero in the backward extrapolation in time. In general Eq. ( 13 1 is nonlinear and the 
fixed point can be found by a Newton-Krylov iteration. However, this requires many additional 
LBM evaluations to construct the Jacobian. Similarly, matrix-free methods like GMRES still 
require many matrix- vector products since the spectrum is unfavorable for fast convergence [31 j - 
In [33] the CR-algorithm is combined with Newton's method by performing local updates at the 
ghost points of the hybrid model to reduce the size of the Jacobian. 



3.3 Discussion of existing lifting operators 



The methods discussed in Sections |3.1| and |3.2| are well known methods to construct a lifting 
operator for LBMs. However, each of these methods has some drawbacks. As noted earlier, a 
drawback of the use of the Chapman-Enskog expansion (Section 3.1 1 is the necessity to construct 
the expressions analytically. Therefore its use is limited to a few examples where the expansion is 
known. However, its computational cost is limited to the calculation of the numerical approxima- 
tion of the derivatives. The Chapman-Enskog expansion then becomes a stencil operator and the 
cost of the application grows linearly with the number of points where lifting is required. 
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Figure 2: Sketch of the first few steps of the Constrained Runs algorithm for the LBM with a con- 
stant backward extrapolation in time. The solid line shows the evolution of {p{x, t), <j)(x, t), £(x, t)} 
along the slow manifold for a given grid point x. For a given po we search for the intersection 
of the plane with the slow manifold. We start iterating with po, the known density, and initial 
guesses <fio and £o f° r the missing moments (v°). After each step of the LBM, the density is reset 
to its initial value po but the moments evolve during the LBM time simulation. This results in 
v k , the k-th iterate of the CR-algorithm. This algorithm finds an approximation for the missing 
values (f> and £ on the slow manifold. 



The Constrained Runs scheme (Section 3.2) can be used to approximate these expressions 
numerically However, the lifting method can become computationally expensive since it requires 
many evaluations of the underlying lattice Boltzmann model to construct the Jacobian matrix. 
Even with the matrix-free methods and local updates discussed at the end of Section 3.2 it still 
remains computationally expensive to use in practice, especially in higher dimensional problems. 

As an advantage of the Constrained Runs algorithm, we should note that the lifting error 
can be smaller than the modeling error, the difference in density between the LBM and its PDE 
approximation, by using the CR-algorithm in the hybrid model discussed in Section 2.3 [29 . 



Section pA\ contains a comparison of the computational cost of these existing methods in the sense 
of hybrid models. 

The focus of this paper is to obtain an alternative lifting operator that reduces the computa- 
tional cost but holds the advantage of achieving the modeling error. 



4 Numerical Chapman-Enskog expansion 

In this section, we construct a lifting operator that alleviates the computational expense of the 
CR-algorithm. It combines the ideas of Constrained Runs and the Chapman-Enskog expansion. 
Instead of using Constrained Runs to find for each grid point the missing moments <p and £ of the 
distribution functions, we use Constrained Runs to find the unknown coefficients of the Chapman- 
Enskog expansion. This has several advantages that will be discussed at the end of Section |4.5[ 

The derivations in this section are again based on one-dimensional problems but can easily be 
generalized to more dimensions. 
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4.1 Distribution functions as a series of the density 

This section shows that the solution fi(x, t), i £ { — 1, 0, 1} of a LBM with an infinite domain and 
parameters Ax, At and uj can be written as a series of p{x,t), the macroscopic density. We ini- 
tially characterize distribution functions /, as smooth functions that are sufficiently differentiablc 
functions in time and space which implies that the same holds for the density, a sum of these 
distribution functions. The smoothness condition will be specified below. This condition can be 
justified when the lattice spacing Ax is much bigger than the mean free path [161 [7]. Then the 
distribution functions can be written as 

ft +\ ten *\, d P , a Q2 P , x ° 3 P , d4 P , 

dp > d 2 p d 2 p ,, JN 

+ ^4 + ^W + - + r,i dxm + - (14) 



where 





(15) 



are fixed constants that only depend on uj, Aa; and At. The derivation of this expansion is outlined 
in the remaining of this section. 

Since the functions fi(x,t) are infinitely differentiable, a Taylor expansion can be constructed. 
The distribution functions in point x + iAx, i £ { — 1, 0, 1} at time t + At are given by 

df- d 2 f i 2 Ax 2 df- d 2 f-At 2 d 2 f- 

M x + iAx, t + At) = fi (x, t) + J±iAx + ^ — + At + _ + ^AxAt + 

Combined with the assumption that /j is a solution of the LBE ^ on an infinite domain, we end 
up with 

^ lAxdh i 2 Ax 2 d 2 fi Atdfr At 2 d 2 f t iAxAtd 2 f % 



fi{x,t) = ft q (x,t) 



uj dx 2uj dx 2 uj dt 2uj dt 2 uj dxdt 

With the notation Li for the functional 

iAx d i 2 Ax 2 d 2 At d At 2 d 2 iAxAt d 2 
1 uj dx 2uj dx 2 uj dt 2uj dt 2 uj dxdt ' 

we can rewrite the LBE into a set of three coupled PDEs for the distribution functions 

(l-£ i )f i (x,t)=f° g (x,t), V^e {-1,0,1}, (17) 

that holds for x £] — oo, oof and t £ [0, oof. 

The solution can be found by performing a Picard or fixed point iteration 

= A/W + f-\ (18) 

with initial guess f\ X ' = that results in 

/f = f?=: 9 i, 

/| 2) = 9i + £i9i + &(£i9i), 



/r = 



fe=0 
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with f!> the n-th iterate. This iteration converges if the error between subsequent iterations goes 
to zero. 

In contrast to traditional iterations, which require convergence for any initial guess, Eq. ( 18 1 is 
a fixed point iteration with initial guess zero and a smooth right hand side. It is only necessary to 
show convergence for this particular case. To discuss this convergence we introduce the 2-norm, 
||.||, to show what happens between subsequent iterations. The absolute difference of subsequent 
iterations is given by 

ll/f +1) -/i n) ll = IIA n+ Vl|. 

This goes to zero if /j is smooth enough, implying smoothness on p and gi = f^ q {x,t) such that 
lim„_ i . 00 = 0. This smoothness condition depends on the parameters of the LBM, Ax, 

At and Li, and the derivatives of gi. For example, when gi can be described by a polynomial, we 
have that there exists a k such that for all n > k applies that ||£" +1 gj|| = 0. 

We end up with the series (14) that consists of the vectors of constants given in (151, the 
density and its derivatives. Once the constants are determined, the lifting operator — that is 
necessary to initialize the LBM and to determine the ghost points in the hybrid model — can be 
constructed. How these constants are found is discussed in Sections 14.21 and 14.31 Section l!~2l deals 



with the analytical derivation while Section 4.3 is concerned with the numerical procedure. As a 
surplus, it allows us to find the corresponding macroscopic PDE as outlined in Section [4. 4| 



4.2 Derivation of a lifting operator 



With the help of expansion (14) it is possible to build a lifting operator that constructs the 



distribution functions for a given density. The focus of this section lies in the determination of the 



vectors of constants ( 15 ), the coefficients of such a lifting operator (14). To simplify the discussion 



and notation we limit ourselves to a truncated series 

dp 



(19) 



where a, /3 and 7 are the vectors containing the constants. The method is easily generalized to 
include higher order terms which will be considered in Section [4. 3.2| 
Using the fact that Eq. (Il9l) is valid for every possible grid point 



every possible grid point, we can consider three grid 

points Xj , Xk and xi and set up a linear system for the nine unknowns, namely three vectors each 

containing three constants. Where j, k and I are certain indices determined in a later stage of the 
paper. 



ft 

ft 

0-1 



fo(&k 1 t) - fa q {^k 1 t) 
f-A^k, t) ~ !'_\(x*,t) 

fi(,x,,t) - /«(*,,*) 



(20) 



For a given fi(x,t) where i G { — 1,0,1}, the linear system (20) will give the coefficients a, j3 
and 7. However, linear system ( p0| only delivers the correct coefficients if fi is smooth enough 
such that f1 q satisfies the smoothness condition. This is the case when fi lies on the slow manifold. 
The Constrained Runs algorithm offers a way to reach the slow manifold in an iterative way. 

We combine the ideas of the CR-algorithm to reach the slow manifold and the Chapman- 
Enskog expansion to find the unknown constants on this slow manifold. The numerical procedure 



to do so is given in Section 4.3 
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If a PDE in closed form exists that describes the evolution of p in the form of p t + ap x = Dp x 



then the linear system (20) will be singular. Indeed, the PDE will give a relation between p t , 
p x and p xx in each of the grid points Xj, Xk and xi. As a result, every element in the last three 



columns of the linear system (20) can be written as a linear combination of the first six columns. 



In practice, however, the PDE is only an approximation and the system will be close to singular. 



This clarifies why we do not solve the linear system in ( 20 1 and first focus on one that is not 



close to singular in Section |4~3) Section [4~4] explains why such a PDE exists. 



4.3 Numerical procedure to construct the lifting operator 

From the previous discussion it is clear that the coefficients of the lifting operator can be extracted 
from a linear system once / approaches the slow manifold. Next, the extraction of the coefficients 
is combined with the CR-algorithm to bring / close to the slow manifold. 

4.3.1 Coefficients of the Chapman-Enskog expansion as a fixed point 

This discussion is limited, as in Section [4. 2| to the first few terms of the expansion. The singular 
system can be avoided by taking a series that only contains spatial derivatives. Such a series can 
represent the same state since the time derivative dtp is often related to the spatial derivatives 
through a macroscopic PDE. For example, suppose that a PDE of the form p t + ap x = Dp xx 
describes the behavior of p. It is then possible to eliminate dtp from the expansion. The coefficients 
are then a — a — 7a and (3 — (3 + 7D. The distribution functions are now series with only spatial 
derivatives. 

f( x ,t) = r(x,t) + a^ + ^. (21) 



Remark 1 Rewrite a and /3 as a and j3 but bear in mind that it considers different coefficients 



in Eq. (191 and Eq. (21) 



Again, once the distribution functions are close to the slow manifold, we can extract the 
coefficients a and f3 from the linear system 



/ 



dp{xj) 






d 2 p(x 3 ) 






dx 


dp{xj) 
dx 


dp(xj) 
dx 


dx 2 


d 2 p{xj) 
dx 2 


d 2 p(x 3 ) 
dx 2 


dp(x k ) 






d' 2 p{x k ) 






dx 


dp(x k ) 
dx 


dp(x k ) 
dx 


dx 2 


d 2 p{x k ) 
dx 2 


d 2 p(x k ) 
dx 2 



/ "1 \ 

V £-1 



fi(xj,t) - 
fo{x 3 ,t) - 



n\ Xjl t) 

fo \ x j'V 



t) 



fl(Xk,t) - 

/o(zfc,*) - 
V f-i(xk,t) 



n q {x kl t) 

fo q (xk,t) 

-n(x k ,t) ) 



(22) 



To reach the slow manifold we combine Constrained Runs with the extraction of the coefficients. 
Consider a numerical function h(a, (3; p, m) as described in Function [T] This function takes as 
input a and f3 and as parameters a fixed density p and an integer m, the order of the smoothness 
condition. It first constructs, with this input, a state / with the help of series ( [2l] ). This state is 
then used to perform multiple LBM steps. For each of these steps we can find the corresponding 
moments <j) and £. On the moments, we can use the CR-algorithm (Section 3.2) to find new 
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moments that are closer to the slow manifold by considering the finite difference approximations 
of the m-th order smoothness condition, 

d m+1 6 d m+1 £ 

a —^=Q and 5— f=0. (23) 

dt m+l dt rn+1 ^ ' 

These new moments result in new coefficients a and /3, by applying the linear system Eq. ( |22[ ) on 
the distribution functions fi, corresponding to the new (f> and £ and the given p. 

The idea is now to determine a and j3 such that they are invariant under this numerical function 
h(a, /3; p, m). Indeed, if the initial and final state can be described by the same a and /3 then the 
lifted / is close to the slow manifold since it is a fixed point of the underlying CR-iteration. 

Instead of performing a regular fixed point iteration with h(a, /?; p,m), a Newton iteration 
is used that finds a :— (a, ft) G K 6 such that a = h(a; p,m). This reduces the computational 
cost significantly because the size of the Jacobian system with a and ft is much smaller then the 
Jacobian of the original Constrained Runs algorithm. The latter involves the moments in every 
grid point and this becomes very large. 

Function 1 h(a, (3; p, m) 

Require: Guess on coefficients a, j3, given density p, order m to use in Eq. (23). 
l: Construct lifting operator / = f e i + a§£ + (Eq. Q). 

2: Compute corresponding moments (f> and £ by applying Eq. (J3|. 

3: Perform m + 1 LBM time steps to compute d di . m +f — and d dtm +i — by using forward finite 
difference formulas. This results in new moments <fi and £. {find moments closer to slow 
manifold} 

4: Revert back to distribution functions / by applying Eq. ([3|. 
5: Select grid points Xj, Xk and xi to construct linear system ( p2| . 
6: Solve the system for a and /?. 
7: return a, j3. 

The numerical function h(a, /3; p, m) has a density p(x, 0) as a parameter and the solution for 
the coefficients is independent of its choice of p. The coefficients are determined by functional 



Li (161 that only depends on constants u), the spatial grid size Ax and time step At. Since the 



coefficients do not depend on time, we can choose an arbitrary density p(x, 0) such that Eq. (22 1 
is easily solvable. 



The choice of the grid points Xj and Xk in (22 1 should be such that the condition number of 



the matrix is optimal. In addition, the spatial derivatives that are considered needs to exist and 
should not become zero during the LBM evolution since otherwise we would end up with singular 
linear systems. 

Furthermore, the test domain used in the LBM inside the function h(a, /3; p,m) can be sig- 
nificantly smaller than the domain of the original LBM problem. A smaller test domain will not 
affect the constant coefficients of the lifting operator. However, it should use the same Ax and At 
as the LBM of interest since the coefficients depend on the chosen spacings in space and time. The 
choice for the test domain, density and indices is further discussed in Section [5. 1| for the considered 
model problem. 

4.3.2 Higher order versions 

There are two ways to increase the accuracy. First, more terms in the expansion can be considered 
such that more derivatives of the density are taken into account. Second, we can enforce a higher 
order smoothness in the CR-algorithm. Both methods are outlined below. 

The proposed method can easily be extended by considering more terms with higher order 



derivatives in the truncated series (21). For example, consider the expansion 
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that now requires the determination of more coefficients that are found by considering — in 
addition to Xj and Xk — additional grid points x\ and x m . This leads to a larger system of 
unknowns but will give better results. 

Higher order smoothness can be enforced on the moments <p and £ as in the CR-algorithm by 
considering a higher order m in Eq. ( 23 1 . This requires more LBM steps and uses a higher order 
finite difference formula to estimate the derivatives in time. 

For further conclusions and results higher order derivatives and higher order smoothness are 
taken into account. 



4.4 Derivation macroscopic PDE 



Next, we derive from Eq. (14), the macroscopic PDE by summing over the velocities. 
J2i fi( x , t) — p( x t t) = J2i fifo, t) results in a macroscopic PDE for the density. 



Using 





d 2 P 

dxdt 



dx 2 




Series ( 14) derived in this setting leads to the classical Chapman-Enskog expansion 28, 5] that 
we obtained in Section 2.2 Indeed, Eq. (Tl4|) is written as 



£i9i + A 9i 



because of the application of the fixed point iteration. When the series is truncated after the 
second order spatial derivative and the first order time derivative, we end up with 



MM) = \p 



iAx dp 
3w dx 



3w 2 



6u> 



d 2 P _ Atdp 
dx 2 3lu dt 



Summing fi(x,t) over i € { — 1,0, 1} we obtain the same macroscopic diffusion PDE Q. Substi- 
tuting this PDE in the series to remove the time derivative leads to the classical Chapman-Enskog 
expansion in Eq. ([5]). 



This macroscopic PDE was used for the removal of the time derivative in Eq. ( 19 ) and its 
replacement with Eq. (|21|). Furthermore, it is important to note that the macroscopic PDE is not 



necessarily of the reaction-diffusion prototype. Truncating ( 14) after more terms and taking more 



derivatives into account results in a better output but it will lead to the term ^§ which gives a 
less comfortable macroscopic PDE. 

This relation to the macroscopic PDE can now be integrated in the numerical Chapman-Enskog 
method. Once the fixed point described in Section [4.3. 1| is found, we have a and /3 that lifts p 
to the distribution functions close to the slow manifold. By performing two more LBM steps, ^ 
can be calculated by using a forward finite difference formula. System ( 20 ) can be applied to find 



the vectors of constants of this larger system that include the time derivative. There are now two 
possibilities: either the resulting system is non-singular and it can be solved for the coefficients 
and only an approximate PDE can be found as is considered above. Or it is too singular to be 
solved accurately but then the PDE can be extracted from the nullspace of the system. 



Let us first discuss the situation where the matrix in (20) is non-singular. The system can 



then be solved for a, (3 and 7. The approximate PDE can be determined by summing over the 
obtained coefficients. 

= J2t a * dp _ J2i & d 2 p 
dt Ei 7; dx J2i 7i dx 2 ' 

This PDE is only approximate. Otherwise, if it would hold exactly, the system would be singular 
as expected. 
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For a singular system, we know that one or more of the eigenvalues will be zero with a cor- 
responding null eigenvector. Focusing on the null eigenvector v = {v\,V2, ■ ■ • ,^9}, we know that 



Av = with A the matrix in system (20 1. Using this, we obtain 



+ + dJ ^Av 7 = 0, 

Ox ox at 

from which we conclude that the resulting PDE looks like 

dp{xj) v\ dp(xj) V4 d 2 p{xj) 
dt Vf dx V7 dx 2 

Remark that same PDE will be found when considering the equation in grid points Xk and x\ 
instead of Xj. 

4.5 Algorithm for lifting operator and macroscopic PDE 

The results of the previous sections are now combined in an algorithm that delivers a lifting 
operator and an approximate macroscopic PDE. This can be used, for example, to construct 
a hybrid model. The pseudocode is presented in Algorithm [2] while the complete algorithm is 
presented below. The algorithm starts by searching for the lifting operator on the basis of the 
spatial derivatives. Thereafter, it inserts time derivatives and calculates the coefficients of the 
macroscopic PDE. 

Start with an initial guess for {a, (3,8, e} in Eq. (24). Apply Function [I] h(a, (3, S, e; p, m) for a 



given p and a certain m for the order of smoothness. This results in coefficients {a,/3, <5, e} that 
represent distribution functions closer to the slow manifold. The lifting operator is constructed at 
this point. 

When these distribution functions are found based on the spatial derivatives only, we still need 
to determine the corresponding PDE by considering the null eigenvector or by a summation of the 
coefficients as discussed in Section |4.4| By performing two extra LBM steps — to estimate the 
time derivative with a forward finite difference formula — the coefficient 7 belonging to the time 
derivative of the expansion below can be numerically calculated. 

f = f eq + a — + f3 — +5— +e— + — 
dx dx 2 dx 3 dx 4 dt 

Since the PDE can be obtained from the numerically constructed distribution functions, the PDE 
obtained through the Chapman-Enskog expansion does not need to be obtained analytically. 

Remark 2 Note that one can also consider f eq (x,t) = np(x,t) and determine the constants of 
vector k in a similar setting by using an extra grid point to obtain a larger system of unknowns. 

Remark 3 In this paper we have chosen to use the same Ax and At in the LBM as in the PDE. 
However, their stability properties may be different. Our specific LBM simulation is stable in the 
2-norm when < ui < 2 J23f . However, it is not necessary that the macroscopic PDE, when it 
is discretized with the same Ax and At and forward Euler, is also stable. Indeed, when uj — > 
for a fixed Ax and At the resulting diffusion coefficient D grows, see Q, and this can lead to an 
instability. 

The pseudocode of the numerical Chapman-Enskog expansion as a lifting operator is given in 
Algorithm [2] together with the determination of the transport coefficients of the PDE to construct 
a hybrid model. The proposed algorithm has several advantages. In contrast to the Chapman- 
Enskog expansion no analytical work is required. Compared to the Constrained Runs algorithm it 
significantly reduces the number of unknowns in the lifting since it only needs to find the coefficients 
(vectors of constants) rather than the full state of the distribution functions. Furthermore, it can 
be done off-line before the calculations. Indeed, once the coefficients are found they can be reused 
every time step to realize the lifting. As an extra surplus, the PDE can be determined to construct 
hybrid models. 
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Algorithm 2 Pseudocode numerical Chapman-Enskog expansion 

Require: Test domain that defines p(x, 0), initial guess a = {a,/3,6,e} = zeros(12,l) and a 
user-defined tolerance tol, parameter m for higher order smoothness, 
repeat 

a k+1 = a k — (j(a fe )) 1 /i(a fe ; p,m) with h defined in Function]!] 
until ||a fc+1 - a k \\ < tol. 

Result for coefficients {a, 0, S, e} belonging to the spatial derivatives. 

Distribution functions / = f eq + + ftgjz + <5g^§ + egpr closer to the slow manifold. 

Perform 2 more LBM steps to determine ^ numerically by a forward finite difference formula. 

Construct system (20) — including higher order spatial derivatives — to achieve coefficients 

{a,/3,<5,e,7}. 

Determine coefficients of PDE by using the nullspace of the system or by summation of the 
coefficients (Section 4.4). 

return Lifting operator / = f eq + aff + /30 + 5^ + e0 and macroscopic PDE. 



5 Numerical Results 

The new lifting operator is now illustrated in several examples. First we benchmark its accuracy 



we 



against a reference solution that is reconstructed. This is done in Section 5.1 In Section 5.2 
recall the one-dimensional hybrid model of Figure [TJ Two-dimensional problems are discussed 
in Section |5.3| The important comparison of the additional required LBM steps to perform the 
lifting in a hybrid model is presented in Section [5. 4| 



5.1 Numerical comparison of different lifting operators 

The proposed lifting operator can be tested against a reference distribution function f c . This 
reference solution is calculated by performing 1000 lattice Boltzmann steps starting from an initial 
state that corresponds to the equilibrium distribution function of a given density p. 

The lifting operator can now be evaluated by restricting the reference distribution function f c 
to its density p = fi{x,t) and lift it back to a distribution function / by using the proposed 
lifting operator. The resulted / will be compared with f c with the help of the 2-norm ||/ — / c ||. 

Example 1 The considered model problem has the following parameters for a one- dimensional 
domain of length L. 

L = 10, n = 200, Aa;=-=0.05, At = 0.001, p(x, 0) = e^"^, ui = 0.9091. 

n 

For these parameters the classical Chapman-Enskog expansion predicts a diffusion coefficient D = 1 
(Eq.m 



To reproduce the numerical results linked to Example [T] we include some extra information on 



how to determine the indices of (22). As mentioned in Section 4.3.1 we can choose an arbitrary 



initial density p(x, 0) and a test domain to determine the constants of the lifting operator ( plj ). 
For example, consider p(x, 0) = x + l/2x 2 for unknowns a, (3 and p(x, 0) = x + l/2x 2 + \/Qx 6 for 
unknowns a, f3, S. x is defined by the spatial nodes of the test domain defined below. Furthermore, 
the test domain reduces the computational expense compared to the actual spatial domain. 

Consider the domain parameters of Example [TJ Since we know that the constants are only 
affected by these space and time steps, we should consider — together with p(x, 0) — a test 



domain with the same step sizes since the vectors of constants ( 15 ) are affected by these choices. 
The test domain is of length L tes t = 3 such that n tes t — 60 since Ax — 0.05. This number of grid 
points will make it possible to choose the indices Xj , Xk, ... such that system ( 22 ) is not close to 



singular. We can return now to the question which indices should be used in system (22). Focus 
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on the fact that we do not want an effect of wrongly chosen boundary conditions in the smaller 
test domain. The grid points should be taken far enough from the edges and in points such that 
the system (22) does not become singular. The indices can be, for example, 10, 20, 30, . . . spread 
over the test domain of 60 grid points. 

Note that one can also focus on local updates around the considered grid points Xj, Xk, 
When the number of iterations needed in Newton's method are known, one knows how many LBM 
steps will be performed to find the new coefficients a and j3. Then the size of the test domain can 
be shrunk to a smaller domain around Xj, Xk, ■ ■ ■■ This is the same idea as used in |33| to perform 
local updates for the CR-algorithm. 

To compare the proposed lifting operator with the existing ones discussed in Section [3] the 
results for ||/ — f c \\ of the different lifting operators are included in Tables [l] [2] and [3] 

Table [l] contains results obtained with the analytical Chapman-Enskog expansion as a lifting 
operator. The first column gives the order of the expansion and the second column shows ||/— / c ||. 
As expected a better accuracy is obtained with higher order expansions. For example, by taking 
the third derivative of the density into account an error of 1.24e-5 is achieved. 



Table 1: The error ||/ — f c \\ is presented to test the exact Chapman-Enskog expansion as a lifting 
operator. The reference solution is obtained after 1000 LBM steps in the model problem from 
Example [I] 



Construction lifting operator: 
exact Chapman-Enskog expansion 


11/ -/ell 


f = f eq 

/ = / e « + cwct§£ 

f f eg j_ _ dp _i_ a dp 
J — J ~T~ inexact q x \ Pexact q x 2 

p req j_ _ dp , r> dp 2 , c 9 3 p 
J — J i Pexact q x ~T~ Pexact g~£Z "exact 


0.0388 
5.2341e-004 
2.7570e-005 
1.2439e-005 



In Table [2] we use the Constrained Runs algorithm of various orders of accuracy to numerically 
lift the density to distribution functions. Different types of backward extrapolation are listed in 
the first column of the table while the corresponding 2-norm ||/ — f c \\ is described in the second 
column. There we see that very accurate results can be found for the higher order versions. Note 
that these methods find for each grid point the moments (f> and £ of the distribution functions. 
Together with p, the corresponding distribution functions are found by Eq. ([3]). Since this gives 
a local solution it can give accurate results. The last column contains ||/ — f c \\ when some extra 
advection effect is included, which shows similar results as the pure diffusion problem. 



Table 2: The error ||/ — f c \\ with the Constrained Runs algorithm (combined with Newton's 
method) for various orders of accuracy as a lifting operator. The reference solution is obtained 
for the model problem in Example [T] by performing 1000 LBM time steps before restricting and 
lifting. The last column contains results when an extra advection effect of a = 0.66 is included — 
which changes the used equilibrium distribution functions as noted in Eq. (Hm. 



Construction lifting operator: 
extrapolation CR-algorithm 


H/-/c|| 

pure diffusion D = 1 


11/ -/c|| 

plus advection a = 0.66 


Constant 


0.0010 


0.0014 


Linear 


1.3578e-006 


1.7927e-006 


Quadratic 


2.9359e-009 


3.9069e-009 


Cubic 


9.0125c-012 


1.1898c-011 



Table [3] shows the results with the proposed numerical Chapman-Enskog expansion as a lifting 
operator. We clearly see that taking more terms in the expansion, i.e. more derivatives of p(x, i) 
in the lifting, leads to a better lifting operator. In the same table we show the results with higher 
order smoothness conditions by using Eq. (23) with higher order m. As in the CR-algorithm, 
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higher order smoothness does result in a significant improvement. The accuracy increases to 
9.45e-ll when up to the sixth spatial derivative is taken into account. Including advection in this 
table will also show similar results but these are not added. 



Table 3: The error ||/ — / c || with the numerical Chapman-Enskog expansion as a lifting opera- 
tor. The reference distribution function f c is obtained by performing 1000 LBM time steps with 
parameters listed in Example [T] Each of the lifting operators is calculated as a fixed point for the 
coefficients. The first table shows the results with one LBM step before updating the moments, 
implying an update of the coefficients. There are results listed for increasing number of terms in 
the expansion, implying an increasing number of considered coefficients. The second table shows 
the same results where two LBM steps are used to estimate the smoothness. The third and fourth 
table show the results with a quadratic, respectively cubic computation of the finite difference 



approximation in Eq. ( 23 ) . The constant coefficients found via the numerical procedure are com- 



pared to those found exactly with the classical Chapman-Enskog expansion in columns 3, 4 and 
5. 



Construction lifting operator: 
Numerical Chapman-Enskog 
constant computation of fixed point 
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OX 
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ox 
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3.7153e-016 
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linear computation of fixed point 
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1.2647e-016 
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1.1455e-014 
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1.1928e-015 
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5.2341e-004 
2.7570e-005 
6.8677e-007 
4.1047e-008 
1.4995e-009 
9.4492e-011 



7.9768e-016 
4.1708e-015 
1.2825e-014 
3.9984e-014 
6.9980e-013 
2.2682e-011 



/ 

1.0779e-014 
3.5932e-014 
7.9750e-014 
1.6787e-012 
5.5927e-011 



/ 
/ 

1.0803e-005 
1.0803e-005 
1.0803e-005 
1.0803e-005 



A comparison of Tables [T] [2] and [3] shows that the proposed numerical lifting operator leads 
to better results than the analytically found Chapman-Enskog expansion. As can be seen is the 
Constrained Runs algorithm a good lifting method, but, as will be discussed in Section |5.4[ the 
computational expense of this method brings down the beauty of it. Table|4]of Section 5.4 contains 
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a comparison of the number of additional LBM steps required for each of the lifting operators. In 
the CR-algorithm this additional cost can be attributed to the construction of the Jacobian matrix. 
These additional LBM steps make the method computationally very expensive. In two dimensions 
this method becomes prohibitive. This makes the numerical Chapman-Enskog lifting operator a 
good alternative to the CR-algorithm that gives a similar accuracy at a limited computational 
cost. 



5.2 One-dimensional test problem 

To compare the results of the numerical Chapman-Enskog expansion with the earlier proposed 
lifting operators discussed in Section [3j Figures [3] and [4] show the results of the absolute difference 
l/\ybrid — Plbm I f° r the exact Chapman-Enskog expansion and those obtained with the Constrained 
Runs algorithm. p hybrid is the density of the hybrid model and /9 LBM the density of a full LBM. 
Plbm 1S t ne reference solution to compare the hybrid solution with. It considers a LBM on the 
whole spatial domain [a, b] with the parameters outlined in Example[l]and the domain represented 
in Figure [T] The lifting operators are used both to initialize the LBM and to find the ghost points 
of the LBM domain. Figure [3] shows the absolute differences by using as lifting operator the 
exact Chapman-Enskog expansion respectively up to zeroth (top left), first (top right), second 
(bottom left) and third order (bottom right). Figure [4] shows the absolute differences with the 
lifting operator based on the Constrained Runs algorithm in combination with Newton's method 
for respectively a constant (top left), linear (top right), quadratic (bottom left) and cubic (bottom 
right) extrapolation in time. These results were obtained in |33j by considering local updates at 
the ghost points of the LBM domain. 

When the numerical Chapman-Enskog expansion (up to the sixth spatial derivative) is used in 
our one-dimensional hybrid model problem, Figure [5] is obtained. Here, we have two possibilities. 
First, act as if we know the PDE (4| obtained from the exact Chapman-Enskog expansion. |p hybrid — 
Plbm I i s gi ven i n the left Figure [5 for which the hybrid domain is shown in Figure [l] and the PDE is 
the one given in Eq. Q. Second, use the PDE that is obtained from the proposed lifting operator 
through summing the proposed lifting operator or considering the nullspace as explained in Section 
With this PDE, the result for \p hybIid — p LBM \ is shown in the right Figure [5] 
As can be seen in Figure [5] a change in the PDE — by considering the PDE obtained through 
the numerical Chapman-Enskog expansion — results in an even smaller modeling error compared 
to the one obtained via the classical Chapman-Enskog expansion. 

Changing the parameters of the model such that advection is included, is considered below. 
The figures show similar results with advection-term a = 0.66. The model problem remains the 
one from Example [I] The only difference is the change in the equilibrium distribution as shown 
in ([6]). Figure [6] contains the comparison results obtained through the CR-algorithm. Figure [7] 
(left) shows |p hybrid — PlbmI with the numerical Chapman-Enskog expansion as a lifting operator 
and the PDE obtained through the Chapman-Enskog expansion while Figure [7] (right) uses the 
PDE obtained from the proposed lifting operator. 



4.4 



5.3 Two-dimensional test problem 

This section generalizes the previous one. Two spatial dimensions are considered. Two-dimensional 
problems can take different discrete sets of velocities into account. In Section [5.3.1| results for D2Q5 
are presented while Section [5. 3. 2| contains results for D2Q9. 



5.3.1 D2Q5 

The hybrid test domain for this section is represented in Figure [8] for D2Q5. Again, the domain 
is split into subdomains. One part of the domain is described by the LBM while another part is 
described by a macroscopic PDE. Example [2] describes the parameters for the model problem in 
this two-dimensional setting. 
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j : Xj,j 6 {0, ... , 199} j : Xj,j 6 {0. . . . , 199} 



Figure 3: \p bybtid — P LB mI after 200 time steps. The difference is also shown at earlier time slots, 
but shifted down for clarity. The lines represent time steps between one and 200. The top 
line corresponds to time step 191 while the bottom line represents time step 11. The lines in 
between correspond to jumps with 10 time steps from 11 to 21, . . .,181 to 191. The domain that 
is considered, is shown in Figure [T] The lifting operator is fi(x,t) — 1/3 p(x,t) (top left), first 
order (top right), second order (bottom left) and third order Chapman-Enskog (bottom right) 
respectively. The lifting operator is used both to find the ghost points of the LBM domain and 
for the creation of the initial state for the LBM region. The model problem parameters are listed 
in Example [l] 
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50 100 150 200 50 100 150 200 

j : xj,j 6 {0, ... , 199} j : Xj,j 6 {0, . . . , 199} 



Figure 4: |p hybrid — P LB m I after 200 time steps where the lifting operator based on the CR-algorithm 
is used in combination with the method of Newton. We show results for constant (top left), linear 
(top right), quadratic (bottom left) and cubic (bottom right) backward extrapolation, respectively. 
The difference is also shown at earlier time slots, but shifted down for clarity. The domain that is 
considered, is shown in Figure [T] The model problem parameters are listed in Example [T] 
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50 100 150 200 50 100 150 200 

j : xj,j 6 {0, ... , 199} j : Xj,j 6 {0,.. ., 199} 

Figure 5: \p hybrid — p LBM \ after 200 time steps in the model problem of Example [I] The domain 
that is considered, is shown in Figure [T] To deal with the initial error and the error in the ghost 
points of the LBM domain the numerical Chapman-Enskog expansion (order spatial expansion 6) 
is used. The considered PDE in the hybrid domain is the analytically known PDE Q in the left 
figure and the one that is obtained from the numerical Chapman-Enskog expansion in the right 
figure. 
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50 100 150 200 50 100 150 200 

j : xj.j 6 {0, ... , 199} j : Xj,j 6 {0. . . . , 199} 



Figure 6: |p hybrid — p LBM I after 200 time steps where the lifting operator based on the CR-algorithm 
is used in combination with the method of Newton. We show results for constant (top left), linear 
(top right), quadratic (bottom left) and cubic (bottom right) backward extrapolation, respectively. 
The difference is also shown at earlier time slots, but shifted down for clarity. The domain that 
is considered, is shown in Figure [l] Model parameters are listed in Example [l] with an extra 
advection coefficient a — 0.66. 




50 100 150 200 50 100 150 200 

j : e {0, . . . , 199} j : x jt j 6 {0, ... , 199} 



Figure 7: \p hybrid — P LB mI a ft er 200 time steps in the model problem of Example [T| with advection 
effect a — 0.66. The domain that is considered, is shown in Figure [T] To deal witn the initial error 
and the error in the ghost points of the LBM domain the numerical Chapman-Enskog expansion 
(order spatial expansion 4) is used. The considered PDE in the hybrid domain is the analytically 
known PDE (7} in the left fi gure and the one that is obtained from the numerical Chapman-Enskog 
expansion in the right figure. 
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Example 2 The considered model problem has the following parameters for a two-dimensional 
domain — described by 5 possible velocity directions (D2Q5) — of length L x L (with n 2 the 
number of grid points). 



L = 10, n = 200, Ax = Ay 



L 



At = 0.0001, 



1.6129. 



For these parameters the classical Chapman- Enskog expansion predicts a diffusion coefficient D = 1 
(Eg. R 



The comparison of \p hyhIid — p LBM I i s represented in Figure [9] for Example [2j The different lifting 
operators are used to obtain distribution functions from a given density. The used lifting operators 
are the equilibrium distribution function in the top left figure, the first order Chapman-Enskog 
expansion in the top right, the second order Chapman-Enskog expansion in the middle left, the 
numerical Chapman-Enskog expansion of spatial order expansion 4 in the middle right and the 
bottom — depending on the used PDE in the hybrid model. 



i e {0,1,..., 4} for D2Q5 

j€{0,...,p} 
k€ {0,...,n-l} 

PDE domain 

Je{p + 1 

ke {0,...,n-l} 
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Figure 8: The two-dimensional spatial domain [a, b] x [a, b] C M 2 in the hybrid model is split into 
[a, l[x [a, b] on which we solve the LBM and [I, b] x [a, b] on which we solve the PDE model. The 
solid points (•) represent the grid for the density p of the discrete PDE, the circles (o) represent 
the LBM variables fi{x, y, t), i € {0, . . . , 4} for D2Q5. The periodic boundaries and the coupling- 
are implemented with ghostcells which are drawn by dashed circles. The density in the ghostcells 
of the PDE domain, in {x p ,yu), k £ {0, ...,n— 1} and (x n ,yk), are found by taking Ylifi m 
(x p ,yj) and (xo,yk), respectively. However, the ghostcells for the LBM domain, in (x—i,yk) and 
(x p+ i, yk), require a lifting operator that lifts p to the distribution functions in these points. 



5.3.2 D2Q9 

This section takes more directions for the velocities into account. Example [3] contains the model 
problem parameters for D2Q9. 

Example 3 The considered model problem has the following parameters for a two-dimensional 
domain — described by 9 possible velocity directions (D2Q9) — of length L x L (with n 2 the 
number of grid points). 

£ = 10, n = 200, Ax = Ay=-, At = 0.00001, u = 1.9531. 
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Figure 9: Comparison between different lifting operators for the model problem presented in Ex- 
ample^ The absolute difference |p hybrid — p LBM | is plotted at time step 200. To deal with the initial 
error and the error in the ghost points of the LBM we use: the equilibrium distribution function 
(top left), the first order Chapman-Enskog expansion (top right), the second order Chapman- 
Enskog expansion (middle left), the numerical Chapman-Enskog expansion (order expansion 4) 
where the PDE in the hybrid domain is the analytically known PDE given in ^ (middle right) 
and the numerical Chapman-Enskog expansion where the considered PDE in the hybrid domain 
is the one that is obtained from the numerical Chapman-Enskog expansion (bottom). 



23 



For these parameters the classical Chapman- Enskog expansion predicts a diffusion coefficient D = 1 

(Eq. 

First consider no advection in the equilibrium distribution functions (a = (0;0)). The results 
for |/3 hybrid — jO LBM | are presented in Figure 10 The used lifting operators are the equilibrium 



distribution (top left), the first order Chapman-Enskog expansion (top right), the second order 
Chapman-Enskog expansion (middle left), the numerical Chapman-Enskog expansion (order ex- 
pansion 4) where the PDE in the hybrid domain is the analytically known PDE given in (|9| 
(middle right) and the numerical Chapman-Enskog expansion where the considered PDE in the 
hybrid domain is the one that is obtained from the numerical Chapman-Enskog expansion (bot- 
tom). 



When advection (a = (1; 0.5)) is included, we end up with Figure 11 for the absolute difference 
IPhybrid - Plbm I when the numerical Chapman-Enskog expansion is used to lift density to distribution 
functions. 



5.4 Analysis of the computational cost of lifting 

This section compares the computational cost of the lifting operators in the one-dimensional test 



problem (Section 5.2) 



The motivation for this paper is to bring down this cost. Especially the Constrained Runs 
algorithm requires many additional LBM steps to lift the density in the ghost points of the LBM 
domain. While numerical Chapman-Enskog only requires a single calculation with a fixed cost 
that can be done off-line before the simulation. This significantly reduces the cost of the lifting. 

A detailed analysis of the lifting cost in terms of additional LBM steps is listed in Table |4j The 
table is an extension of the results of [33 with results for the classical Chapman-Enskog expansion, 
Constrained Runs algorithm combined with Newton's method and the numerical Chapman-Enskog 
expansion. 

It can be seen that the total number of LBM steps for the CR-algorithm are listed per ghost 
point and per time step. The number for the numerical Chapman-Enskog expansion is the total 
for the entire domain and at all time steps since the calculations for the coefficients are done 
off-line. 

In two-dimensional problems the numerical Chapman-Enskog expansion still has a limited 
computational cost. Only a few additional coefficients need to be determined associated with the 
extra spatial derivatives. 

Note that the computational cost of applying the numerical Chapman-Enskog lifting operator 
is the same as applying the analytical Chapman-Enskog operator. For each grid point we need the 
derivatives of p, which can be calculated by finite differences using the densities at neighboring 
grid points. 



6 Conclusions 

This article proposes a numerical lifting operator for lattice Boltzmann models (LBMs) that maps 
a given density to the corresponding distribution functions. This new lifting operator is based on 
the Chapman-Enskog expansion that writes the missing distribution functions as analytical series 
of the density and its derivatives. The coefficients of this expansion are now determined through 
a numerical method, in contrast to the original expansion where they are found analytically. The 
numerical method is based on the Constrained Runs algorithm that relies on the attraction of the 
dynamics toward the slow manifold. 

A systematic numerical comparison of the accuracy and the computational cost between the 
analytical Chapman-Enskog expansion, the Constrained Runs algorithm and the new lifting oper- 
ator is performed in this article. The cheapest way to lift is with the Chapman-Enskog expansion. 
However, the analytical expressions are not always available for the system of interest. An al- 
ternative numerical lifting operator is Constrained Runs (CR), but its computational cost grows 
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Figure 10: |jO hybrid — /3 LBM | after 200 time steps with model problem presented in Example p] 
without advection terms. To deal with the initial error and the error in the ghost points of the 
LBM the equilibrium distribution is used (top left), the first order Chapman-Enskog expansion 
(top right), the second order Chapman-Enskog expansion (middle left), the numerical Chapman- 
Enskog expansion (order expansion 4) where the PDE in the hybrid domain is the analytically 
known PDE given in ([8| (middle right) and the numerical Chapman-Enskog expansion where the 
considered PDE in the hybrid domain is the one that is obtained from the numerical Chapman- 
Enskog expansion (bottom). 
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Figure 11: \p h brid — p LBM | after 200 time steps with model problem presented in Example^ with 
advection a = (1; 0.5). To deal with the initial error and the error in the ghost points of thcLBM 
the numerical Chapman-Enskog expansion (order expansion 4) where the PDE in the hybrid 
domain is the analytically known PDE given in |9| (left) and the numerical Chapman-Enskog 
expansion where the considered PDE in the hybrid domain is the one that is obtained from the 
numerical Chapman-Enskog expansion (right). 



Table 4: Analysis of the computational cost of the lifting operators in terms of the additional 
LBM steps. Both Constrained Runs (CR) and the numerical Chapman-Enskog (NCE) require 
additional LBM steps to lift p to the distribution function /. The table shows these additional 
steps to construct the Jacobian operator for the Newton iteration. The listed values for the 
CR-algorithm even are per ghost point and per time step. While those of the proposed lifting 
operator are the total number for the entire domain and for all time steps together since the 
vectors of constants needs to be determined only once and can be reused throughout the rest of 
the simulation. 



Lifting operator 


Number of 


LBM steps 


Total number 




iterations 


to perform 


of LBM steps 






one iteration 




Exact Chapman-Enskog 


/ 


/ 





CR-algorithm 






per ghost point 


type of extrapolation in time 






per time step 


Constant 


3 


19 x 1 


57 


Linear 


3 


31 x 2 


186 


Quadratic 


3 


43 x 3 


387 


Cubic 


3 


55 x 4 


660 


Numerical Chapman-Enskog 






for entire domain 


with 18 unknowns 






and all time steps 


{a, 13, 8, e, 6,1,} 








Constant 


3/(4) 


19 


57+2=59 


Linear 


3/(4) 


19 x 2 


114+2=116 


Quadratic 


3/(4) 


19 x 3 


171+2=173 


Cubic 


3/(4) 


19 x 4 


228+2=230 
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significantly with the order of accuracy. It needs many additional LBM steps to find the missing 
distribution functions. 

The new result and the main focus of this paper is a numerical lifting method that combines 
the ideas of Constrained Runs and the Chapman-Enskog expansion. Instead of using Constrained 
Runs to find for each grid point the missing moments, we use Constrained Runs to find the 
unknown coefficients of the Chapman-Enskog expansion. This numerical lifting method has several 
advantages. First, it significantly reduces the number of unknowns in the lifting step: we only need 
to find the coefficients rather than the full state f(x,v,t). And secondly, it can be done off-line 
before the calculations. Indeed, once the coefficients are found they can be reused every time step 
and every grid point to realize the lifting, at no significant additional cost. A third advantage is 
that the expansion gives, as a spin-off, the transport coefficients of the macroscopic PDE. 

The new lifting operator, the numerical Chapman-Enskog expansion, is then used in a hybrid 
domain that spatially couples a macroscopic partial differential equation (PDE) with a lattice 
Boltzmann model. This creates a missing data problem at the interfaces since the PDE model has 
too few variables to provide the LBM with the correct boundary conditions. 

The numerical Chapman-Enskog expansion deals with this mismatch in variables. It maps 
the variables of the PDE model to those of the LBM. We evaluate and compare various lifting 
operators. In particular, we have focused on a simple LBM and PDE model discretized with equal 
grid and time steps such that the error created by the coupling can be highlighted. The paper 
presents numerical results both for ID and 2D hybrid domains where part of the LBM domain 
is replaced by the macroscopic PDE. In both cases the error associated with the coupling can be 
made smaller than the modeling error, related to the PDE approximation of the LBM. 

This paper reports on our initial efforts where we have focused on a simple model problem with 
several limiting assumptions. In the model we have assumed an equilibrium distribution function 
that depends only on the local density, while in general it also depends on the local momentum and 
temperature. This limitation can be easily alleviated by considering a Chapman-Enskog expansion 
with a more general equilibrium function. 

A further assumption is that we used the same time and space grid for the PDE and the 
LBM. This choice was made to highlight the error made by the coupling mechanism, the ease of 
implementation and to eliminate the error due to the different discretizations. However, there is 
no reason to prohibit different grid and time spacings. Extra care is then needed to interpolate 
between time and grid spacings. In practice, the grid of the PDE can be further coarsened, 
depending on local discretization errors. Ideally, the hybrid model is embedded in an adaptive 
mesh refinement simulation, where at the finest level a LBM is used. 

We have also kept the boundary between the LBM and the PDE domain fixed during the 
simulation at an arbitrary position. In the future, this boundary should be moved adaptively 
using an accuracy requirement based on the properties of the lifting operator. 

For the model problem with periodic boundary conditions studied in this paper, the Chapman- 
Enskog expansion exists everywhere and we could in principle put the boundary between the PDE 
and the Boltzmann model at any location, provided that we lift accurately. For general Boltzmann 
models, with complicated collision integral operators, such a Chapman-Enskog expansion might 
not exist everywhere in the domain. Then a hybrid model can be constructed where a PDE can 
replace the Boltzmann model only in the regions where the Chapman-Enskog expansion is known 
to exist. 

This situation appears in the modelling of laser ablation where a laser heats a surface that 
consequently melts and evaporates. The escaping plasma plume can be described by a Boltzmann 
equation. Close to the melting surface a complicated non-equilibrium situation appears where 
escaping particles evaporate but particles that impinge on the melted surface condensate. There 
is no Chapman-Enskog expansion that can describe this situation close to the surface. Only away 
from the surface the plasma reaches an equilibrium situation. A hybrid model will then use a full 
Boltzmann model near melt while a reduced PDE model can be used away from the surface. 
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